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Abstract 

'^ ' We introduce the classical stellar atmosphere problem and describe in detail its 

QQ ■ numerical solution. The problem consists of the solution of the radiation transfer 
equations under the constraints of hydrostatic, radiative and statistical equilibrium 

, _ (non-LTE). We outline the basic idea of the Accelerated Lambda Iteration (ALI) 

^-^ • technique and statistical methods which finally allow the construction of non-LTE 

cn . model atmospheres considering the influence of millions of metal absorption lines. 

_ Some applications of the new models are presented. 

o 
(^ 



j-< , 1 Introduction 

O 

c/3 , The quantitative analysis of stellar spectra is one of the most important tools 

^ ' of modern astrophysics. Basically all our knowledge about structure and evo- 

lution of stars, and hence about galactic evolution in general, rests on the in- 
^ ' terpretation of their electromagnetic spectrum. The formation of the observed 

j^ ■ spectrum is usually confined to a very thin layer on top of the stellar core, 

the atmosphere. Spectral analysis is performed by modeling the temperature 
and pressure stratification of the atmosphere and computing synthetic spec- 
tra which are then compared to observation. Fitting synthetic spectra from 
a grid of models yields the basic photospheric parameters, effective tempera- 
ture, surface gravity, and chemical composition. Comparison with theoretical 
evolutionary calculations allows the derivation of stellar parameters like mass, 
radius and total luminosity. 

The so-called classical stellar atmosphere problem considers the transfer of 
electromagnetic radiation, released by interior energy sources, through the 
outermost layers of a star into free space by making three specific physical 
assumptions. At first it is assumed that the atmosphere is in hydrostatic equi- 
librium, thus, the matter which interacts with photons is at rest. Second, the 
transfer of energy through the atmosphere is entirely due to photons, i.e. heat 
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conduction and large scale convection are regarded as negligible (so-called 
radiative equilibrium). The effectiveness of photon transfer depends on the 
total opacity and emissivity of the matter which are strongly state and fre- 
quency dependent quantities. They depend in detail on the occupation density 
of atomic levels which in turn are determined by the local temperature and 
electron density as well as by the radiation field, whose nature is non-local 
in character. The occupation of any atomic level is balanced by radiative and 
collisional population and de-population processes (statistical equilibrium; our 
third assumption), i.e. the interaction of atoms with other particles and pho- 
tons. Mathematically, the whole problem consists of the solution of the ra- 
diation transfer equations simultaneously with the equations for hydrostatic 
and radiative equilibrium, together with the statistical equilibrium, or, rate 
equations. 

A stellar atmosphere is radiating into the circumstellar space and thus ev- 
idently is an open thermodynamic system, hence it cannot be in thermody- 
namic equilibrium (TE) and thus we cannot simply assign a temperature. The 
"Local Thermodynamic Equilibrium" (LTE) is a working hypothesis which as- 
sumes TE not for the atmosphere as a whole but for small volume elements. As 
a consequence, the atomic population numbers are depending only on the lo- 
cal (electron) temperature and electron density via the Saha-Boltzmann equa- 
tions. Computing models by replacing the Saha-Boltzmann equations by the 
rate equations are called non-LTE (or NLTE) models. This designation is un- 
fortunate because still, the velocity distribution of particles is assumed to be 
Maxwellian, i.e. we can still define a local temperature. NLTE calculations 
are tremendously more costly than LTE calculations, however, it is hard to 
predict if NLTE effects are important in a specific problem. Generally, NLTE 
effects are large at high temperatures and low densities, which implies intense 
radiation fields hence frequent radiative processes and less frequent particle 
collisions which tend to enforce LTE conditions. 

Relaxing the LTE assumption leads to the classical model atmosphere prob- 
lem, i.e. solution of the radiation transfer equations assuming hydrostatic, 
radiative and statistical equilibrium. Such models are applicable to the vast 
majority of stars. The numerical problem going from LTE to realistic NLTE 
models has only recently been solved and is the topic of this paper. We now 
have the tools in hand to consider non-classical models, which consider the 
radiation transfer in more general environments, for example in expanding 
stellar atmospheres. This is the topic of another paper in this volume [1] . 

Stellar atmosphere modeling has made significant progress within the recent 
years. This is based on the development of new numerical techniques for model 
construction as well as on the fact that reliable atomic data have become avail- 
able for many species. Of course these achievements go along with a strong 
increase of computing power. Model atmospheres assuming LTE have been 



highly refined by the inclusion of many more atomic and molecular opacity 
sources, however, elaborated numerical techniques for LTE model computa- 
tion are available for many years. The progress is most remarkable in the 
field of NLTE model atmospheres. The replacement of the Saha-Boltzmann 
equations (LTE) by the atomic rate equations (NLTE) requires a different 
numerical solution technique, otherwise metal opacities cannot be accounted 
for at all. Such techniques were developed with big success during the last 
decade, triggered by important papers by Cannon [2] and Scharmer [3]. The 
Accelerated Lambda Iteration (ALT) is the basis of this development. Com- 
bined with statistical methods we are finally able to compute so-called metal 
line blanketed NLTE models (considering many millions of spectral lines) with 
a very high level of sophistication. 

In this paper we discuss the basic ideas behind the new numerical methods 
for NLTE modeling. At first we state the classical model atmosphere problem 
and describe the ALI solution technique. We then focus on the NLTE metal 
line blanketing problem and its solution by the introduction of the superlevel 
concept and statistical methods to treat the opacities (Opacity Sampling and 
Opacity Distribution Functions). Finally we demonstrate successful applica- 
tions of the new models by presenting a few exemplary case studies. 



2 Statement of the problem and overview of the solution method 



In the following text we outline the general stellar atmosphere problem, but 
will discuss various details of numerical implementation as applied to our com- 
puter program PR02. We assume plane parallel geometry, which is well justified 
for most stars because the atmospheres are thin compared to the stellar ra- 
dius. The only parameters which characterize uniquely such an atmosphere are 
the effective temperature {Tes), which is a measure for the amount of energy 
transported through the atmosphere per unit area and time (see Eq. 15), the 
surface gravity (g), and the chemical composition. Generalization to spherical 
symmetry to account for extended (static) atmospheres mainly affects the ra- 
diation transfer equation and is straightforward [4, p. 250]. To construct model 
atmospheres we have developed our program which solves simultaneously a set 
of equations that is highly coupled and non-linear. Because of the coupling, no 
equation is determining uniquely a single quantity - all equations determine 
a number of state parameters. However, each of them is usually thought of as 
determining a particular quantity. These equations are: 

• The radiation transfer equations which are solved for the (angular) mean 
intensities Ji,i = 1, . . . ,NF, on a pre-chosen frequency grid comprising 
NF points. The formal solution is given by J = AS", where S is the source 
function as defined later (Eq. 11). Although A is written as an operator, one 



may think of A as a process of obtaining the mean intensity from the source 
function. 

• The hydrostatic equihbrium equation which determines the total particle 
density A^. 

• The radiative equilibrium equation from which the temperature T follows. 

• The particle conservation equation, determining the electron density Ue- 

• The statistical equilibrium equations which are solved for the population 
densities n^, i = 1, . . . , NL of the atomic levels allowed to depart from LTE 
(NLTE levels). 

• The definition equation for a fictitious massive particle density Uh which is 
introduced for a convenient representation of the solution procedure. 

This set of equations has to be solved at each point rf of a grid comprising 
ND depth points. Thus we are looking for solution vectors 

il)'a = {ni,...,nNL,ne,T,nH,N,Ji,...,JNF), d = l,...,ND. (1) 



The Complete Linearization (CL) method [5] solves this set by linearizing the 
equations with respect to all variables. The basic advantage of the ALI (or 
"operator splitting") method is that it allows to eliminate at the outset the 
explicit occurrence of the mean intensities Jj from the solution scheme by 
expressing these variables by the current, yet to be determined, occupation 
densities and temperature. This is accomplished by an iteration procedure 
which may be written as (suppressing indices indicating depth and frequency 
dependency of variables) : 

J" = A*^" + (A - A*)5"-^ (2) 

This means that the actual mean intensity at any iteration step n is computed 
by applying an Approximate Lambda Operator (ALO) A* on the actual (ther- 
mal) source function S*" plus a correction term that is computed from quanti- 
ties known from the previous iteration step. This correction term includes the 
exact lambda operator A which guarantees the exact solution of the radiation 
transfer problem in the limit of convergence: J = AS. The use of A in Eq. 2 
only indicates that a formal solution of the transfer equation is performed but 
in fact the operator is usually not constructed explicitly. Instead a Feautrier 
solution scheme [4, p. 156] or any other standard method can be employed to 
solve the transfer equation that is set up as a differential equation. 

The resulting set of equations for the reduced solution vectors 

^l^a = {ni,...,nNL,ne,T,nH,N), d=l,...,ND (3) 



is of course still non-linear. The solution is obtained by linearization and it- 
eration which is performed either with a usual Newton-Raphson iteration or 
by other, much faster methods like the quasi-Newton or Kantorovich vari- 
ants [6,7]. The first model atmosphere calculations with the ALI method were 
performed by Werner [8]. 

Another advantage of the ALI method is that the explicit depth coupling of 
the solution vectors Eq. 1 through the transfer equation can be avoided if one 
restricts to diagonal (i.e. local) approximate A-operators. Then the solution 
vectors Eq. 3 are independent from each other and the solution procedure 
within one iteration step of Eq. 2 is much more straightforward. Depth cou- 
pling is provided by the correction term that involves the exact solution of the 
transfer equation. The hydrostatic equation which also gives an explicit depth 
coupling, may be taken out of the set of equations and can - as experience 
shows - be solved in between two iteration steps of Eq. 2. Then full advantage 
of a local ALO can be taken. 

The linearized system may be written as 

V', = V'^ + S^f^a (4) 

where Vd i^ ^^e current estimate for the solution vector at depth d and 5i/'d 
is the correction vector to be computed. Using a tri-diagonal A* operator 
the resulting system for 6il)^ is - like in the classical CL scheme - of block 
tri-diagonal form coupling each depth point d to its nearest neighbors d±l: 

Id^'^d^i + (^d^'^d + "d^V'd+i = Crf. (5) 

The quantities ck, /3, 7 are {NNxNN) matrices where NN is the total number 
of physical variables, i.e., A^A^ = NL + 4, and 0^ is the residual error in the 
equations. The solution is obtained by the Feautrier scheme. With starting 
values Di = /3^^(-q;i) and vi = /3j~^Ci we sweep from the outer boundary of 
the atmosphere inside and calculate at each depth: 



^d={f3d + ld^d-iy\cd-Ja^d-i)- (6) 

At the inner boundary we have Dnd = and sweeping back outside we 
calculate the correction vectors, first ^i/jjvd = ^nd and then successively 
6ipd = TydS'4^d+i + Vrf. As already mentioned, the system Eq. 5 breaks into 
ND independent equations (5i/j^ = /3^ c^ {d = 1, . . . , ND) when a local A* 
operator is used. The additional numerical effort to set up the subdiagonal 
matrices and matrix multiplications in the tri-diagonal case is outweighed by 



the faster global convergence of the ALI cycle, accomplished by the explicit 
depth coupling in the linearization procedure [9]. 

The principal advantage of the ALI over the CL method becomes clear at this 
point. Each matrix inversion in Eq. 6 requires (A^L + 4)^ operations whereas in 
the CL method {NL + NF + 4)^ operations are needed. Since the number of 
frequency points NF is much larger than the number of levels NL, the matrix 
inversion in the CL approach is dominated by NF . 

Recent developments concern the problem that the total number of atomic 
levels tractable in NLTE with the ALI method described so far is restricted to 
the order of 250, from our experience with PR02. This limit is a consequence 
of the non-linearity of the equations, and in order to overcome it, measures 
must be taken in order to achieve a linear system whose numerical solution is 
much more stable. Such a pre-conditioning procedure has been first applied 
in the ALI context by Werner & Husfeld [10]. More advanced work achieves 
linearity by replacing the A operator with the \E' operator (and by judiciously 
considering some populations as "old" and some as "new" ones within an ALI 
step) which is formally defined by writing 

J^ = <^uVy , i-e. ^^ = K/Xu , (7) 

where the total opacity Xv (as defined in Sect. 3.7) is calculated from the 
previous ALI cycle. The advantage is that the emissivity rj^ (Sect. 3.7) is linear 
in the populations, whereas the source function S^ is not. Hence the new 
operator \E' gives the solution of the transfer problem by acting on a linear 
function. This idea is based on Rybicki & Hummer [11] who applied it to the 
line formation problem, i.e. restricting the set of equations to the transfer and 
rate equations and regarding the atmospheric structure as fixed. Hauschildt et 
al.[l,12] generalized it to solve the full model atmosphere problem. In addition, 
splitting the set of statistical equations and solving it separately for each 
chemical element means that now many hundreds of levels per species are 
tractable in NLTE. A very robust method and fast variant of the ALI method, 
the ALI/CL hybrid scheme, allows for the linearization of the radiation field 
for selected frequencies [13], but it is not implemented in PR02. 



3 Basic equations 

3.1 Radiation transfer 



Any numerical method requires a formal solution (i.e. atmospheric structure 
already given) of the radiation transfer problem. The radiation transfer at any 



particular depth point can be described by the following equation, formally 
written for positive and negative fi (which is the cosine of the angle between 
direction of propagation and outward directed normal to the surface) sepa- 
rately, i.e. for inward and for outward directional intensities / with frequency 
u: 

±f^^^^^ = S,-I,{±fi), fie [0,1]. (8) 

T^ is the optical depth (which can be defined via the column mass m that is 
used in the other structural equations and later introduced in Sect. 3.3.2 by 
dTi, = dmXu/p, with the mass density p) and S^, is the local source function. 
Introducing the Feautrier variable 

u,^ = {lM + I,{-fi))/2 (9) 

we obtain the second-order form [4, p. 151]: 

f) If 

^^'^ = ^^,-Su, f^e[0,i]. (10) 



We may separate the Thomson emissivity term (scattering from free electrons, 
assumed coherent, with cross-section o"e) from the source function so that 

Su = Sl + HeCTeJu/Xu, (H) 



where Si is the ratio of thermal emissivity to total opacity as described in 
detail below (Sect. 3.7): S'^ = rju/xu- Since the mean intensity is the angular 
integral over the Feautrier intensity the transfer equation becomes 

yu r. 2 = "i^M - ^u / "i^M dp- (12) 



Thomson scattering complicates the situation by the explicit angle coupling 
but the solution can be obtained with the standard Feautrier scheme. Assum- 
ing complete frequency redistribution in spectral lines [4, p. 29], no explicit 
frequency coupling occurs so that the parallel solution for all frequencies en- 
ables a very efficient vectorization on the computer. 

The following boundary conditions are used for the transfer equation. At the 
inner boundary where the optical depth is at maximum, r = Tmaxj we have 

~ K^l ~ "^unijia&x) (13) 




where we specify /j^ = /zy(+Ai, Tmax) from the diffusion approximation: 

1+ =B + ^^ - iU) 



Bjj is the Planck function and Ti the nominal (frequency integrated) Eddington 
flux: 

n = anT^jATi (15) 



with the Stefan-Boltzmann constant a^. At the outer boundary we take Ty = 
Tmin = iTT'iXu/'^P, assumlug that X is a linear function of m for m < mi. Since 
Tmin 7^ 0, it is uot cxactly valid to assume no incident radiation at the stellar 
surface. Instead we specify I^^ = I^^—ji^ Tmin) after Scharmer & Nordlund [14]: 

Cm = 'S',,(rmm)[l - exp(-rmm//i)] (16) 



which follows from Eq. 8 assuming S{t) = S{t^u^) for r < Tmin- Then we get 

f/^^T^I = Uy,.{rnim) - I'f,- (17) 



' ' nim 



The boundary conditions are discretized performing Taylor expansions which 
yield second-order accuracy [4, p. 155]. 



3.2 Statistical equilibrium 



The statistical equilibrium equations are set up according to [4, p. 127]. The 
number of atomic levels, ionization stages and chemical species, as well as 
all radiative and collisional transitions are taken from the input model atom 
supplied by the user (Sect. 6.3). Ionization into excited states of the next 
ionization stage is allowed for. Dielectronic recombination and autoionization 
processes can also be included in the model atom. 



3.2.1 Rate equations 

As usual the atomic energy levels are ordered sequentially by increasing exci- 
tation energy, starting with the lowest ionization stage. Then for each atomic 



level i of any ionization stage of any species the rate equation describes the 
equilibrium of rates into and rates out of this level: 






The rate coefficients Pij have radiative and collisional components: Pij ■ 
Rij + Cij. Radiative upward and downward rates are respectively given by: 



oo 



(TijiU] 



R^. =47r / -^j^J^du (19) 





Photon cross-sections are denoted by o"jj(z/). {ni/rij)* is the Boltzmann LTE 
population ratio in the case of line transitions: Qi/gjexpl—huij/kT), where 
the Qij are the statistical weights. The LTE population number of a particular 
level is defined relative to the ground state of the next ion, so that in the case 
of recombination from a ground state nf we have by definition {ni/rij)* = 
ne(f)i{T) with the Saha-Boltzmann factor 

(PAT) = 2.07 ■ lO-^(i3^J.-3/2^huJkT (21) 

9i 

where hi^i is the ionization potential of the level i. Care must be taken in 
the case of recombination from an excited level into the next low ion. Then 

{ni/rijY = ne(pi ■ (pt/^f^j- 

Dielectronic recombination is included following [15]. Assuming now that j is 
a ground state of ion k, then the recombination rate into level i of ion k — 1 
via an autoionization level c (with ionization potential huc, having a negative 
value when lying above the ionization limit) is: 

R, - ^«.*/.^--"'"-^? (l + ^) . (22) 

The reverse process, the autoionization rate, is given by: 

4n^e^ 1 

Rij = -, ficJ- (23) 

hmc Vc 

The oscillator strength for the stabilizing transition (i.e. transition i— >c) is 
denoted by /jc, and J is the mean intensity averaged over the line profile. 



The program simply takes J^ from the continuum frequency point closest 
to the transition frequency, which is reasonable because the autoionization 
line profiles are extremely broad. The population of autoionization levels is 
assumed to be in LTE and therefore such levels do not appear explicitly in 
the rate equations. 

The computation of collisional rates is generally dependent on the specific ion 
or even transition. Several options, covering the most important cases, may 
be chosen by the user. 



3.2.2 Abundance definition equation 

The rate equation for the highest level of a given chemical species is redundant. 
It is replaced by the abundance definition equation. This equation simply 
relates the total population of all levels of a particular species to the total 
population of all hydrogen levels. Summation over all levels usually includes 
not only NLTE levels but also levels which are treated in LTE, according to 
the specification in the model atom. Denoting the number of ionization stages 
of species k with NION{k), the number of NLTE and LTE levels per ion with 
NL{1) and LTE{1), respectively, we can write: 



NION{k) 

E 
1=1 



'NL(l) 



E ^kli + 



j=l 



LTE{1) 

E 



n 



kli 





'NL{H) LTE{H) 


= yk 


^ ni+ y^ n\ + np 




i=l j=l 



■ (24) 



On the right hand side we sum up all hydrogen level populations including the 
proton density Up, and Hk is the number abundance ratio of species k relative 
to hydrogen. 



3.2.3 Charge conservation 

We close the system of statistical equilibrium equations by invoking charge 
conservation. We denote the total number of chemical species with NATOM, 
the charge of ion / with q{l) (in units of the electron charge) and write: 



ATOA 


f NION{k) 


'NL{1) LTE{1) 


>; 


>; 1(1) 


E ^kli+ E ^kli 


fc=i 


1=1 


i=l i=l 



rie 



(25) 
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3.2-4 Complete statistical equilibrium equations 

We introduce a vector comprising the occupation numbers of all NLTE levels, 
n = (ni, . . . ^riNLj. Then the statistical equilibrium equation is written as: 

An = b. (26) 



The gross structure of the rate matrix A is of block matrix form, because tran- 
sitions between levels occur within one ionization stage or to the ground state 
of the next ion. The structure is complicated by ionizations into excited levels 
and by the abundance definition and charge conservation equations which give 
additional non-zero elements in the corresponding lines of A. 



3.3 Radiative equilibrium 



Radiative equilibrium denotes the fact that the energy transport is exclusively 
performed by photons. It can be enforced by adjusting the temperature strati- 
fication either during the linearization procedure or in between ALI iterations. 
In the former case a linear combination of two different formulations is used 
and in the latter case a classical temperature correction procedure (Unsold- 
Lucy), generalized to NLTE problems, is utilized. The latter is particularly 
interesting, because it allows to exploit the blocked form of the rate coefficient 
matrix. This will enable an economic block-by-block solution followed by a 
subsequent Unsold-Lucy temperature correction step. On the other side, how- 
ever, this correction procedure may decelerate the global convergence behavior 
of the ALI iteration. 



3.3.1 Differential and integral forms for linearization procedure 

The two forms of writing down the radiative equilibrium condition follow from 
the postulation that the energy emitted by a volume element per unit time is 
equal to the absorbed energy per unit time (integral form): 



xAS,-J,)du = Q, (27) 





where scattering terms in x ^"^^ Sy cancel out. This formulation is equiva- 
lent to invoking flux constancy throughout the atmosphere (differential form) 



11 



involving the nominal flux Ti (Eq. 15): 

fj-ifMdiy-n = (28) 

^''^ 

where fu is the variable Eddington factor, defined as 

1 1 

U = I^^Uy^ djij I Uy^, djj, (29) 



and computed from the Feautrier variable u^^ (Eq. 9) after the formal solution. 
As discussed e.g. in [16] the differential form is more accurate at large depths, 
while the integral form behaves numerically better at small depths. Instead 
of arbitrarily selecting that depth in the atmosphere where we switch from 
one formulation to the other, we use a linear combination of both constraint 
equations which guarantees a smooth transition with depth, based on physical 
grounds [7,17]. Before adding up both equations we have to take two measures. 
At first we divide Eq. 27 by the absorption mean of the opacity, kj, for scahng 
reasons: 

oo oo 

Kj = — / KyJudv with J = J,ydiy, (30) 





where Hu = Xu — neCTe is the true opacity without electron scattering. Then 
we multiply Eq. 28 with a similar average of the diagonal elements of the A* 
matrix: 

oo 

A*j = - f A*J^du. (31) 





These two steps determine the relative weight of both equations in a particular 
depth. Numerical experience shows that it is necessary to damp overcorrec- 
tions by adding the following term, which is computed from quantities of the 
previous iteration step and which vanishes in the limit of convergence, to the 
right hand side of Eq. 28: 




F, = [A^r I ^^iUu) du - A^jH (1 - A}) 
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(32) 
last iterate 



We write the equation of radiative equilibrium in its final form: 

— fxAS.-Ju)diy + A''jfj—{f,J,)diy-A'}n-Fo = 0. (33) 



We note that explicit depth coupling is introduced by the differential form 
Eq. 28 through the derivative d/dry even if a purely local A* operator is 
used. Therefore the linearization procedure can no longer be performed in- 
dependently at each depth point and the question becomes relevant at which 
boundary to start with. Numerical experience shows that it is essential to 
start at the outer boundary and to continue going inwards. If a tri-diagonal 
operator is used, nearest neighbor depth coupling is introduced anyhow. The 
program user can choose either the linear combination Eq. 33 or the purely 
integral form Eq. 27, the latter may be necessary to start the iteration under 
certain circumstances. The linear combination, however, is found to give a 
much faster convergence behavior. 



3.3.2 Unsold- Lucy temperature correction procedure 

Closely following Lucy [18] (but avoiding the Eddington approximation and 
using variable Eddington factors instead) and generalizing to NLTE one can 
derive for each depth point a temperature correction AT to be applied to the 
actual temperature in order to achieve flux constancy. Using dx = —dr^/xui 
the zeroth momentum (i.e. angle averaged form) of the radiation transfer Eq. 8 
is: 

^ = X. {S, - J.) ■ (34) 

dx 



with Ju from Eq. 2 and the Eddington fiux H^. In the LTE case with elec- 
tron scattering, S^ can be written as the sum of a thermal and a scattering 
contribution: 

S, = ^B, + ^^ J,. (35) 

Xu Xv 



In the NLTE case we formally write in analogy: 

S, = '^B, + ^ J, (36) 

Xy Xy 
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with quantities k^ and 'ju which can be freely evaluated but which are not 
independent of each other, since •jt, must be expressed by 

Iv = 7 (37) 



in order to yield S,y on the r.h.s. of Eq. 37. With this substitution Eq. 34 reads: 

-^ = K^B, - ix. - 7.) J.- (38) 



Integrating over frequencies, the condition of flux conservation then reads: 

^ = ^J-B^O (39) 



where we used the following definitions for k,, np, and dr: 

OO 1 °° 

K = - J iXu - lu)Ju du = - J {KuiJu - S^) + K^B^)du, (40) 



OO 

— / K^Bj^dv and dr = —Rpdx. (41) 

D J 



np ^ 



Since we can choose k,^ freely, we can define which opacities shall contribute, 
finally resulting in a favorable scaling of factors in Eq. 43. Usually we start 
with all processes included in k^ to begin with moderate corrections. Fol- 
lowing Hauschildt (priv. comm.) one can optionally exclude bound-bound or 
bound-free transitions which is necessary if strong lines or continua dominate 
numerically the radiative equilibrium in optically thin regions. Note that this 
measure does not affect the solution in the case of convergence, but only the 
convergence rate. Without such an acceleration, the Unsold-Lucy procedure 
may run into pseudo-convergence. 

Integrating the first momentum of the radiation transfer equation over fre- 
quency we obtain: 

— — = - — Ti. with Rh = — / f^uHi, dv. (42) 

dr Kp H J 





Using Eq. 39 and the depth integrated form of Eq. 42 we proceed as described 
by Lucy [18]. We finally obtain, with frequency averaged Eddington factors / 
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and h as well as ks defined in analogy to Eq. 30, the temperature correction 
at any depth: 



AT ^ 



AanT^ Kp 



.J J - .sS + J ^ -AH dm + ^(Q^ j 



(43) 



where A7i is the difference between the actual and the nominal Eddington 
flux. In practice it is useful to accelerate this procedure by extrapolating the 
last, say, ten corrections. 

The Unsold-Lucy procedure provides model atmospheres with a relative devi- 
ation from the flux constancy smaller than 10^^ which is a factor of ten better 
when compared to the procedure employing Eq. 28. Due to the decoupling of 
the temperature from the statistical equilibrium the Unsold-Lucy procedure 
is numerically much more stable allowing to calculate models which other- 
wise failed to converge. The price is a slower overall convergence of the ALI 
iteration by a factor of two. 



3.4 Hydrostatic equilibrium 

We write the equation for hydrostatic equilibrium as [4, p. 170]: 

where g is the surface gravity and m the column mass. P is the total pressure 
comprising gas, radiation and turbulent pressures, so that: 




An 7 1 \ 

NkT + — j UJ, dv + -pt;Lb = 9 (45) 



with Boltzmann's constant k and the turbulent velocity fturb- The hydrostatic 
equation may either be solved simultaneously with all other equations or sep- 
arately in between iterations. The overall convergence behavior is usually the 
same in both cases. If taken into the linearization scheme and a local A* opera- 
tor is used then, like in the case of the radiative equilibrium equation, explicit 
depth coupling enters via the depth derivative d/dm. Again, solution of the 
linearized equations has to proceed inwards starting at the outer boundary. 
The starting value in the first depth point (subscript rf = 1) is 



CO 



2 \ c I p, ) 
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where h^ is the variable Eddington factor denoting the ratio of H^/J,, at the 
surface, kept fixed during hnearization. 



3.5 Particle conservation 



The total particle density A^ is the sum of electron density plus the population 
density of all atomic states, LTE and NLTE levels. We may write down the 
particle conservation equation in the following form that contains explicitly 
only the hydrogen population numbers: 



N = ne + 



'NL(H) LTE(H) 

J2 ni+ Y. < + % 



NATOM 

E y^- (47) 

fc=i 



3. 6 Fictitious massive particle density 



A fictitious massive particle density hh is introduced for notational conve- 
nience. It is defined by 

NATOM NATOM 

nH = {N- Ne) E ^kVk/ E Vk- (48) 

k=l k=l 



The mass of a chemical species in AMU is denoted by ruk- Introducing the 
mass of a hydrogen atom mn, we may simply write for the material density 

p = nnmH- (49) 

3.7 Opacity and emissivity 



Thermal opacity and emissivity are made up by atomic radiative bound- 
bound, bound-free and free-free transitions. For each chemical species we com- 
pute and sum up: 



NION 

1=1 



E E ^l^^lA^) \^l^ - n..-e-^(^--^)/'=^) (50) 

j=i j>i \ 9ij J 

NL{1) NL{1+1) 
i=l j>i 
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^NL{1+1) LTE{1+1) 



where the total opacity includes Thomson scattering, i.e. Xu = ^^u + loco's, and 



Vi^ 



2hu^/c^ 



NION 

^ E 



■AfL(0 AfL(0 

E E ffH^«i('^)"-/i^e 



9li_-h(u-u,,)/kT) 



(51) 



AfL(/) AfL(/+l) 

E E 



+ E E ^H^i^iAi^Xe-'''^''' 



^NL{1+1) LTE{1+1) 



The first index of variables marked with two indices denotes the ionization 
stage and the second one denotes the ionic level. Thus au^i^i^ki^) denotes the 
cross-section for photoionization from level i of ion / into level k of ion / + 1. 
The double summation over the bound-free continua takes into account the 
possibility that a particular level may be ionized into more than one level of 
the next high ion. Again, note the definition of the LTE population number 
n*j in this case, which depends on the level {I + l,k) of the parent ion: 



n^i = ni+i^kne(t>ii 



(t>i^ 



1,1 



n+i,k 



(52) 



Note also, that the concept of LTE levels (whose population densities do enter, 
e.g. the number or charge conservation equations) in the atomic models of 
complex ions is therefore not unambiguous. The present code always assumes 
that LTE levels in the model atoms are populated in LTE with respect to the 
ground state of the upper ion. 

The source function used for the approximate radiation transfer is the ratio 
Tjy/ Ky, thus, excludes Thomson scattering. For the exact formal solution of 
course, the total opacity Xv i^^ the expression Eq. 11 includes the Thomson 
term {rieCJe)- 



3.8 Atomic level dissolution by plasma perturbations 



As high-lying atomic levels are strongly perturbed by other charged particles in 
the plasma they are broadened and finally dissolved. This effect is observable 
by line merging at series limits and has to be accounted for in line profile 
analyses. Moreover, line overlap couples the radiation field in many lines and 
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flux blocking can strongly affect the global atmospheric structure. Numerically, 
we treat the level dissolution in terms of occupation probabilities, which for 
LTE plasmas can be defined as the ratio of the level populations to those in 
absence of perturbations. A phenomenological theory for these quantities was 
given in [19]. The non-trivial generalization to NLTE plasmas was performed 
by Hubeny et al.[20]. In practice an individual occupation probability factor 
(depending on T, Ug, and principal quantum number), is applied to each atomic 
level which describes the probability that the level is dissolved. Furthermore, 
the rate equations Eq. 18 must be generalized in a unique and unambiguous 
manner. For details see [20]. As an example. Fig. 1 shows these occupation 
probabilities for hydrogen and helium levels as a function of depth in a white 
dwarf atmosphere. 



1.0 



'0.5 



0.0 - 




2-4 

log column mass / g cm ' 

Fig. 1. Occupation probabilities of energy levels of ionized helium (left panel) and 
hydrogen (right panel) as a function of depth in a white dwarf atmosphere. Levels 
with high principal quantum number n are already dissolved in the upper atmo- 
sphere (left boundary of panels). At the inner boundary even all low- lying levels are 
essentially dissolved. Atmospheric parameters are Tcfj=100 000K, log (7 =7.5, and 
H/He=0.1%. 



4 The Accelerated Lambda Iteration (ALI) 



In all constraint equations described above the mean intensities J^, are sub- 
stituted by the approximate radiation field Eq. 2 in order to eliminate these 
variables from the solution vector Eq. 1. In principle the approximate lambda 
operator may be of arbitrary form as long as the iteration procedure con- 
verges. In practice however an optimum choice is desired in order to achieve 
convergence with a minimum amount of iteration steps. The history of the 
ALOs is interesting and was summarized in detail by Hubeny [21]. Of utmost 
importance were two papers by Olson and collaborators [22,23] who overcame 
the major drawback of early ALOs, namely the occurrence of free parameters 
controlling the convergence process, and who found the optimum choice of 
ALOs. Our model atmosphere program enables the use of either a diagonal or 
a tri-diagonal ALO, both are set up following [23]. 



4-1 Diagonal (local) lambda operators 



In this case the mean intensity Jd at a particular depth d in the current 
iteration step is computed solely from the local source function S^ and a 
correction term AJ^, the latter involving the source functions (of all depths) 
from the previous iteration. Dropping the iteration count and introducing 
indices denoting depth points we can rewrite Eq. 2: 

Jd = A*ddSd + A.U (53) 



In the discrete form we now think of A* as a matrix acting on a vector whose 
elements comprise the source functions of all depths. Then A^ ^ is the diagonal 
element of the A* matrix corresponding to depth point d. Writing A*j^^ = B^ 
(for numerical computation see Eq. 57 below) we have a purely local expression 
for the mean intensity: 

Jd = BaSd + A Jrf (54) 

4-2 Tridiagonal (non-local) lambda operators 



Much better convergence is obtained if the mean intensity is computed not 
only from the local source function but also from the source function of the 
neighboring depths points. Then the matrix representation of A* is of tri- 
diagonal form and we may write 

Jd = Cd^iSd-i + BdSd + Afi+iSd+i + AJd (55) 



where Cd-i and Ad+i represent the upper and lower subdiagonal elements 
of A* and Sd±i the source functions at the adjacent depths. In analogy the 
correction term becomes 

AJd = Ad^'Sd' - {Cd-iSd-i + BdSd + Ad+iSd+i). (56) 



Again all quantities for the computation of AJd are from the previous iteration 
and the first term denotes the exact formal solution of the transfer equation. 
We emphasize again that the actual source functions in Eq. 55 are computed 
from the actual population densities and temperature which are unknown. 
We therefore have a non-linear set of equations which is solved by either a 
Newton-Raphson iteration or other techniques, resulting in the solution of a 
tri-diagonal linear equation of the form Eq. 5. 
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As was shown in [22] the elements of the optimum A* matrix are given by the 
corresponding elements of the exact A matrix. The diagonal and subdiagonal 
elements are computed from [23]: 











with Arrf_i = (r^i — Td-i) / ji. At large optical depths with increasing Ar steps 
(the depth grid is equidistant in logr) the subdiagonals A^+i and Cd^i vanish 
and the diagonal Bd approaches unity, resembling the fact that the radiation 
field is more and more determined by local properties of the matter. At very 
small optical depths all elements of A* vanish, reflecting the non-localness of 
the radiation field in this case. 



4-3 Acceleration of convergence 



PR02 allows usage of an acceleration scheme to speed up convergence of the 
iteration cycle Eq. 2. We implemented the scheme originally proposed by Ng 
[24,25]. It extrapolates the correction vector 511)^ from the previous three it- 
erations. From our experience the extrapolation often yields over-corrections 
resulting in alternating convergence or even divergence. And usually the ap- 
plication of a tri-diagonal ALO results in a satisfactorily fast convergence so 
that the acceleration scheme is rarely used. 



5 Solution of the non-linear equations by iteration 



The complete set of non-linear equations for a single iteration step Eq. 2 com- 
prises at each depth the equations for statistical, radiative, and hydrostatic 
equilibrium and the particle conservation equation. For the numerical solution 
we introduce discrete depth and frequency grids. The equations are then lin- 
earized and solved by a suitable iterative scheme. Explicit angle dependency of 
the radiation field is not required here and consequently eliminated by the use 
of variable Eddington factors. Angle dependency is only considered in the for- 
mal solution of the transfer equation. The program requires an input model 
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atmosphere structure as a starting approximation together with an atomic 
data file, as well as a frequency grid. Depth and frequency grids are therefore 
set up in advance by separate programs. 



5. 1 Discretization 



A depth grid is set up by an auxiliary program which computes, starting from 
a gray approximation, a LTE continuum model using the Unsold-Lucy temper- 
ature correction procedure. In this program depth points are set equidistantly 
on a logarithmic (Rosseland) optical depth scale. The user may choose the in- 
ner and outer boundary points and the total number of grid points (typically 
90). The converged LTE model (temperature and density structure, given on 
a column mass depth scale) is written to a file that is read by PR02. The NLTE 
code uses the column mass as an independent depth variable. 

The frequency grid is established based upon the atomic data input file (see 
Sect. 6.3). Frequency points are set blue- and redward of each absorption edge 
and for each spectral line. Gaps are filled up by setting continuum points. 
Finally, the quadrature weights are computed. The user may change default 
options for this procedure. Frequency integrals appearing e.g. in Eq. 33 are 
replaced by quadrature sums and differential quotients involving depth deriva- 
tives by difference quotients. 



5.2 Linearization 



All variables x are replaced by x — ^ x + (5x where 5x denotes a small per- 
turbation of X. Terms not linear in these perturbations are neglected. The 
perturbations are expressed by perturbations of the basic variables: 

5x = ^ST + ^5n^ + ^5N + ^5nn + E ^5n,. (58) 

oT oue oN oriH j^i oni 



As an illustrative example we linearize the equation for radiative equilibrium. 
Most other linearized equations may be found in [8]. Assigning two indices {d 
for depth and i for frequency of a grid with NF points) to the variables and 
denoting the quadrature weights with Wi Eq. 33 becomes: 



E Wi{—^[5Sdi - 5Jdi\ + 5xdASdi - Jdi\) 
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NF 



Xdi 



NF 



Wi 



2_.Wi^—{Sdi — Jdi) ~ A*j2. -T {fdi-Jdi — fd-l,iJd~l,i) 



(59) 



j=i 



Note that we do not linearize Atj. Because of this, convergence properties 
may be significantly deteriorated in some cases. Perturbations SSdi, Sxdi are 
expressed by Eq. 58, and the perturbation of the mean intensity Jdi is, accord- 
ing to Eq. 55, given through the perturbations of the source function at the 
actual and the two adjacent depths: 



SJdi — Cd-i,i5Sd^i,i + BdiSSdi + Ad-^-i^iSSd+i,i 



(60) 



where A, B, C are the A matrix elements from Eq. 57. The 5Jd-i,i involve the 
term Cd-2,i5Sd-2,i which is neglected because we only want to account for 
nearest neighbor coupling. We write 5Sdi with the help of Eq. 58 and observe 
that for any variable z 



dS, 



di 



dz 



1 (dridi 
Xdi \ 



dz 



Sdi 



dx. 



di 



dz 



(61) 



Derivatives of opacity and emissivity with respect to temperature, electron 
and population densities are computed from analytical expressions (see e.g. 
[26,27]). We finally get from Eq. 59: 



^rr i^ WidSd^^ 

(>id-i,i \ y. —- — 7^ — Xdi^d-i,i 
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= r.h.s. (62) 

Curly brackets {■ ■ ■} denote terms that are similar to those multiplied with 
the perturbations of the temperature. Instead of partial derivatives in respect 
to T, they contain derivatives in respect to ne and the populations n;. They 
all represent coefficients of the matrices ck, /3, 7 in Eq. 5. 



5.3 Newton-Raphson iteration 



As described in Sect. 2 the linearized equations have a tri-diagonal block- 
matrix form, see Eq. 5. Inversion of the grand matrix (= T sized {NN -ND) x 
(A^A^ ■ ND), i.e. about 10^ x 10^ in typical applications) is performed with a 
block-Gaussian elimination scheme, which means that our iteration of the non- 
linear equations represents a mult i- dimensional Newton-Raphson method. The 
problem is structurally simplified when explicit depth coupling is avoided by 
the use of a local ALO, however, the numerical effort is not much reduced, 
because in both cases the main effort lies with the inversion of matrices sized 
A"A^ X A^A^. The Newton-Raphson iteration involves two numerically expensive 
steps, first setting up the Jacobian (comprising ck, /3, 7) and then inverting it. 
Additionally, the matrix inversions in Eq. 6 limit their size to about A^A^ = 150 
because otherwise numerical accuracy is lost. Two variants recently introduced 
in stellar atmosphere calculations are able to improve both, numerical accuracy 
and, most of all, computational speed. 



5.4 Alternative fast solution techniques for non-linear equations: Broyden- 
and Kantorovich-variants 



Broyden's variant [28] belongs to the family of so-called quasi-Newton methods 
and it was first used in model atmosphere calculations in [6,29,30]. It avoids 
the repeated set-up of the Jacobian by the use of an update formula. On top 
of this, it also gives an update formula for the inverse Jacobian. In the case 
of a local ALO the solution of the linearized system at any depth is 



Si^ = l3^'c. (63) 



Let /3ji. ^ be the k-th iterate of the inverse Jacobian, then an update can be 
found from: 



-1 , (sfc- ^fcVfc)^( s^/3,^^ 



Pk+i = Pk + ^r^-1 , (64) 
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where denotes the dyadic product and where we have defined: 

Sfc = (5-0;. solution vector of preceding hnearization 

y^ = c,fc+i — Ck difference of actual and preceding residuum. 

The convergence rate is super-linear, i.e. slower than the quadratic rate of the 
Newton-Raphson method, but this is more than compensated by the tremen- 
dous speed-up for a single iteration step. It is not always necessary to begin 
the iteration with the calculation of an exact Jacobian and its inversion. Ex- 
perience shows that in an advanced stage of the overall (ALI-) iteration Eq. 2 
(i.e. when corrections become small, of the order 1%) we can start the lin- 
earization cycle Eq. 64 by using the inverse Jacobian from the previous overall 
iteration. Computational speed-up is extreme in this case, however, it requires 
storage of the Jacobians of all depths. 

More difficult is the application to the tri-diagonal ALO case. Here we have 
to update the grand matrix T which, as already mentioned, is of block tri- 
diagonal form. We cannot update their inverse, because it is never computed 
explicitly. Furthermore we need an update formula that preserves the block 
tri-diagonal form which is a prerequisite for its inversion by the Feautrier 
scheme Eq. 6. Such a formula was found by Schubert [31]: 

rr rp , (yfc-TfcSfc)0S^ 

Tk+l=i-k-\ ZJfrz (65) 

where Sfc = Zs^ with the structure matrix Z as defined by: 

1 if Tijj^O 



Zij 



if Tij = 0. 



The vectors s^ and y^ are defined as above but now they span over the quan- 
tities of all instead of a single depth point. With this formula we obtain new 
submatrices q:,/3,7 and c with which the Feautrier scheme Eq. 6 is solved 
again. This procedure saves the computation of derivatives. Another feature 
realized in our program also saves the repeated inversion of q = {fB^j^ + j^^Dd-i) 
by updating its inverse with the Broyden formula Eq. 64. Similar to the diag- 
onal ALO case it is also possible to pass starting matrices from one overall 
iteration Eq. 2 to the next for the update of T and the matrix q~^. In both 
cases the user specifies two threshold values for the maximum relative cor- 
rection in Sil^ which cause the program to switch from Newton-Raphson to 
Broyden stages 1 and 2. During stage 1 each new overall cycle Eq. 2 is started 
with an exact calculation and inversion of all matrices involved and in stage 
2 these matrices are passed through each iteration. 
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Another variant, the Kantorovich method was recently introduced into model 
atmosphere calculations [7]. It is more simple and straightforward to imple- 
ment. This method simply keeps fixed the Jacobian during the linearization 
cycle and it is surprisingly stable. In fact it turns out to be even more stable 
(i.e. it can be utilized in an earlier stage of iteration) than the Broyden method 
in the tri-diagonal ALO case. The user of PR02 may choose this variant in two 
stages in analogy to the Broyden variant. It was found that in the stage 2 it is 
necessary to update the Jacobian, say, every 5 or 10 overall iterations in order 
to prevent divergence. 



6 NLTE metal line blanketing 



Despite the capacity increase for the NLTE treatment of model atmosphere 
problems by introducing the ALI method combined with pre-conditioning 
techniques, the blanketing by millions of lines from the iron group elements 
arising from transitions between some 10^ levels could only be attacked with 
the help of statistical methods. These have been introduced into NLTE model 
atmosphere work by Anderson [32,33]. At the outset, model atoms are con- 
structed by combining many thousand of levels into a relatively small number 
of superlevels which can be treated with ALI (or other) methods. Then, in or- 
der to reduce the computational effort, two approaches were developed which 
vastly decrease the number of frequency points (and hence the number of 
transfer equations to be solved) to describe properly the complex frequency 
dependence of the opacity. These two approaches have their roots in LTE mod- 
eling techniques, where for the same reason statistical methods are applied for 
the opacity treatment: The Opacity Distribution Function (ODF) and Opac- 
ity Sampling (OS) approaches. Both are based on the circumstance that the 
opacity (in the LTE approximation) is a function of two only local thermo- 
dynamic quantities. Roughly speaking, each opacity source can be written in 
terms of a population density and a photon cross-section for the respective 
radiative transition: 

>^u ~ ni(Tiu{v) 

In LTE the population follows from the Saha-Boltzmann equations, hence 
rii = ni{ne,T). The OS and ODF methods use such pre-tabulated (on a very 
fine frequency mesh) Ku{ne-, T) during the model atmosphere calculations. The 
NLTE situation is more complicated, because pre-tabulation of opacities is 
not useful. The population densities at any depth now also depend explic- 
itly on the radiation field (via the rate equations which substitute the TE 
Saha-Boltzmann statistics) and thus on the populations in each other depth 
of the atmosphere. As a consequence, the OS and ODF methods are not ap- 
plied to opacity tabulations, but on tabulations of the photon cross-sections 



25 



(t(z/). These do depend on local quantities only, e.g. line broadening by Stark 
and Doppler effects is calculated from T and rif.. In the NLTE case the cross- 
section takes over the role which the opacity played in the LTE case. So, 
strictly speaking, the designation OS and ODF is not quite correct in the 
NLTE context. 

The strategy in our code is the following. Before any model atmosphere calcu- 
lation is started, the atomic data are prepared by constructing superlevels, and 
the cross-sections for superlines. Then these cross-sections are either sampled 
on a coarse frequency grid or ODFs are constructed. These data are put into 
the model atom which is read by the code. The code does not know if OS or 
ODFs are used, i.e. it is written to be independent of any of these approaches. 



6. 1 Model atoms for iron group elements 



The large number of atomic levels in a single ionization stage is grouped into 
a small number of typically 10-20 superlevels or, energy bands. Grouping is 
performed by inspecting a level diagram (Fig. 2) which shows the number of 
levels (times their statistical weight) per energy bin as a function of excitation 
energy. Gaps and peaks in this distribution are used to define energy bands. 
Each of these bands is then treated as a single NLTE level with suitably 
averaged statistical weight and energy. All individual lines connecting levels 
out of two distinct bands are combined to a band-band transition with a so- 
called complex photon cross- sect ion. This cross-section essentially is a sum of 
all individual line profiles which however conserves the exact location of the 
lines in the frequency spectrum. This co-addition is performed once and for 
all and on a very fine frequency mesh to account for the profile shape of every 
line, before any model atmosphere calculation begins. These complex cross- 
sections (examples are seen in the top panels of Figs. 3 and 4) are tabulated 
and later used to construct ODFs or to perform OS for the model calculations. 

Each of the model bands L is treated as one single NLTE level with an average 
energy E^ and statistical weight G^ which are computed from all the individual 
levels {El, Qi) within a particular band: 

E. = T. Eig: h g: G, = e^^"^' E 9*. (66) 

leL I UL leL 

where g* = a^gie~^'^''^* . T* is a characteristic temperature, pre-chosen and 
fixed throughout the model calculations, and at which the ionization stage in 
question is most strongly populated. Energy levels of all iron group elements s 
in the same ionization stage contribute to these model bands according to their 
abundance a^. All individual line transitions with cross-sections (T,„ between 
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Fig. 2. Energy distribution of statistical weights of the iron group elements in ion- 
ization stage VI. Individual energy levels are grouped into bands (horizontal lines) 
and merged into superlevels with an average energy. 

two model bands L and U are combined to one complex band-band transition 
with a cross-section a^u as described by: 
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(67) 



(f^i'i-'iu — T^) is the normalized profile of an individual line. This means that all 
individual lines are correctly accounted for in a sense that their real position 
within the frequency spectrum is not affected by the introduction of atomic 
model bands. The complex cross-sections (each possibly representing many 
thousand individual lines) are computed in advance of the model atmosphere 
calculations on a fine frequency grid with a resolution smaller than one thermal 
Doppler width (typically 0.1 Az/^). This is done at two values for the electron 
density (0 and 10^^ cm^^) and the NLTE code accounts for depth dependent 
electron coUisional broadening by interpolation. Individual line photon cross- 
sections are represented by Voigt profiles including Stark broadening. 

CoUisional excitation rates between atomic model bands are treated with a 
generalized Van Regemorter [34] formula: 



with 
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(68) 
iu)/kT* (69) 



where P{x) = max[5f,0.276e^ii^i(a;)]. Eh is the ionization potential of hydro- 



27 



gen (in electron volts), Ei is the first exponential integral and g a constant 
depending on the ionic charge. The T^u involve the f-values of all individual 
lines and they are computed together with the radiative cross-sections. Third 
degree polynomials in log kT* are fitted to F^^ and the coefficients are written 
into the atomic input data file for the NLTE code. 

Photoionization cross-sections for iron group elements have numerous strong 
resonances that are difficult to deal with. As a first approximation one can cal- 
culate hydrogen-like cross-sections a,,;,/ for the individual levels / and combine 
them to a complex ionization cross-section for every model band: 

^B.M = e-^-/'^*E^X./M- (70) 

UL 



This cross-section is stored in a file and read by the code. Other data to be used 
alternatively (available e.g. from the Opacity Project) may easily be prepared 
and stored in such a file by the user. For collisional ionization one may select 
Seaton's [35] formula with a mean (hydrogen-like) ionization cross-section. 



6.2 OS and ODF approaches 



The OS or, alternatively, ODF approaches are introduced merely in order to 
save computing time during the model atmosphere calculations. In principle it 
is possible to proceed directly with the complex cross-sections constructed as 
described above. However, this would require a very fine frequency mesh over 
the entire spectrum in order to discretize the cross-sections in a similar detailed 
manner, resulting in some 10^ frequency points. Since computation time scales 
linearly with the number of frequency points in the ALI formalism, a reduction 
to some thousand or ten thousand points easily reduces the computational 
effort by an order of magnitude. 

Opacity Sampling is the more straightforward approach. The fine cross-section 
is sampled by a coarse frequency grid and the resulting coarse cross-section 
is used for the model calculation (Fig. 3). Individual lines are no longer ac- 
counted for in an exact way, but this is not necessary in order to account for 
the line blanketing effects, i.e. effects of metal lines on the global atmospheric 
structure like surface cooling and backwarming of deeper layers. A high res- 
olution synthetic spectrum can be obtained easily after model construction 
by performing a single solution of the radiation transfer equation on a fine 
frequency mesh. 

The quality of the sampling procedure can be checked by a quadrature of the 
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Fig. 3. More than 1000 line transitions between individual levels from two spe- 
cific superlevels (Fig. 2) are co-added to a complex photon cross-section resolved by 
330 000 frequency points (top panel). Sampling of this cross section with 300 points 
results in the cross section shown in the bottom panel. 

cross-section on the frequency grid (with weights Wk)'- 
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(71) 



Renormahzation may be performed if necessary. This reduction of the cross- 
sections by samphng is also performed before the model calculations begin. 

The alternative way is the construction of Opacity Distribution Functions 
(or, more correctly, cross-section distribution functions). Each complex cross- 
section is re-ordered in such a way that the resulting ODF is a monotonous 
function (see Fig. 4, middle panel). The resulting smooth run of the cross- 
section over frequency can be approximated by a simple step function with 
typically one dozen of intervals. This cross-section is then fed into the model 
atmosphere code which can use a coarse frequency mesh to appropriately 
incorporate the ODFs. In order to avoid unrealistic systematic effects, however, 
the cross-section bins within each ODF are re-shuffled randomly (bottom panel 
of Fig. 4). 

Many numerical tests concerning model atom construction with superlevels 
were performed by studying the effects of details in band definition and widths. 
Also, the resulting model atmospheres using ODF and OS approaches were 
compared and generally, good agreement was found [36]. 
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Fig. 4. Construction of a cross-section distribution function and its representation 
by a step function (middle panel) from a portion of a complex cross-section (top). 
In the bottom panel a randomized arrangement of the interval steps is shown. 

6.3 Atomic data and model atoms 



All recent progress in stellar atmosphere modeling would have been impossible 
without the availability of atomic input data. Major data sources were put at 
public disposal by the Opacity Project [37] and the work by Kurucz [38]. These 
sources provide energy levels, transition probabilities, and bound-free photon 
cross-sections. The Iron Project [39] is delivering electron collision strengths 
which are important for NLTE calculations and which were hardly available 
up to now. We cannot over-emphasize these vital contributions to our work. 

The atomic species that are to be included in the model atmosphere calcu- 
lations are entirely determined by an atomic data input file. For each ion- 
ization stage of any chemical element the user defines atomic levels by the 
ionization potential and statistical weight and assigning a name (character 
string) to them. These level names are used to define radiative and collisional 
bound-bound and bound-free transitions among the levels as well as free- 
free opacities. The declaration of such transitions is generally complemented 
by a number which specifies the formula which PR02 shall use to calculate 
cross-sections for the rates and opacities. Depending on the formula chosen 
by the user, additional input data are occasionally expected, such like oscilla- 
tor strengths or line broadening data. Alternatively, photon cross-sections for 
lines and continua may be read from external files whose names need to be 
declared with the definition of the transitions. Construction of model atoms 
involving large datasets, e.g. from the Opacity Project, is automated and re- 
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quires a minimum of work by the user. The interested reader is referred to a 
comprehensive User's Guide for PR02 available from the authors [40]. 



7 Application to hot compact stars 



One important motivation for developing and applying the ALT method for 
stellar atmospheres was the unsolved problem of NLTE metal line blanketing 
in hot stars. We want to focus here on two topics which highlight the successful 
application of the new models. The first concerns the Balmer line problem 
which until recently appeared to be a fundamental drawback of NLTE models. 
The second example describes the abundance determination of iron group 
elements in evolved compact stars by constructing self-consistent models which 
can reproduce simultaneously the observed spectral properties of white dwarfs 
and subdwarf O stars (sdO) from the optical region through the extreme 
ultraviolet regime. 



7.1 Balmer line profiles 



Fitting synthetic profiles to observed Balmer lines is the principal ingredient of 
most spectroscopic analyses. The so-called Balmer line problem represents the 
failure to achieve a consistent fit to the hydrogen Balmer line spectrum of any 
hot sdO star whose effective temperature exceeds about 70 000 K. Results of 
Teff determinations drastically differ, up to a factor of two, depending on which 
particular line is fitted. This problem was uncovered a few years ago during 
NLTE analyses of very hot subdwarfs and central stars of planetary nebulae 
[41]. Since then, it cast severe doubt upon NLTE model atmosphere analysis 
techniques as a whole. With new available models computed with the ALI 
method we were able to demonstrate that the problem is due to the neglect 
or improper inclusion of metal opacities [42]. We showed that the Balmer line 
problem can be solved when surface cooling by photon escape from the Stark 
wings of lines from the C, N, O elements is accounted for (see Figs. 5 and 6). 



7.2 Heavy element abundances in sdO stars and white dwarfs 



The optical spectra of hot white dwarfs and sdO stars are dominated by helium 
and/or hydrogen lines. Metals are highly ionized and their spectral lines are 
almost exclusively located in the UV and extreme UV regions. High resolution 
spectroscopy with the International Ultraviolet Explorer (lUE) has revealed a 
wealth of spectral features from iron and nickel which, however, could not be 



31 



90 



80 



70 



60 



S 50 



40 



30 




H-He model 

plus CNO with Doppler profiles 
plus CNO with Stark profiles 
H-He LTE model 



-5 -4 -3 -2 -1 

log column mass / g cm"^ 

Fig. 5. Temperature structures of three NLTE model atmospheres with increas- 
ing degree of sophistication. A LTE model is shown for comparison (triangles). 
While pure H-He models show a high-temperature plateau at the surface, our most 
sophisticated model shows a monotonous temperature run (full line). Formation 
depth intervals for selected lines are represented by horizontal bars. Teff=82 000K, 
log (7 =6.2, solar abundance ratios. 
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Fig. 6. Line profile fits to the subdwarf O star BD-|-28°4211. Two sets of synthetic 
profiles are plotted, namely that from the pure H-He model structure shown in 
Fig. 5 (dashed) and that from the model with CNO elements and Stark broadened 
profiles. The Balmer line problem is occurring with the former set and essentially 
disappears with the latter. Note in particular that the Ha emission core is perfectly 
matched. 

analyzed because of the lack of appropriate NLTE calculations. First attempts 
for quantitative analyses were performed with line formation calculations on 
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pre-specified temperature and pressure model structures which in turn were 
obtained from simphfied LTE or NLTE models, i.e. disregarding metal line 
blanketing effects [43,44]. Subsequently fully line blanketed LTE models were 
employed, however, NLTE effects turned out to be non-negligible [45-49] . Our 
latest models [50] include 1.5 million lines from the iron group elements, which 
are taken from Kurucz's [38] line list. As an example for the quality of the 
fits we can achieve. Fig. 7 shows a portion of the UV spectrum of the sdO 
star Feige 67 and the best fitting model. The derived abundances suggest that 
radiative levitation is responsible for the extraordinarily high heavy element 
abundances in these stars. 




1310 1315 1320 

A / A 



Fig. 7. Model fit (solid line) to the lUE spectrum of an sdO star. Iron and nickel are 
overabundant by factors of 10 and 70, respectively, relative to solar values. Vertical 
bars in the lower half of the panel indicate spectral line positions of four different 
Fe and Ni ions and the bar heights correspond to the relative log gf values (from 

[36]). 



8 Conclusion 



We have described in detail the numerical solution of the classical model at- 
mosphere problem. The construction of metal line blanketed models in hy- 
drostatic and radiative equilibrium under NLTE conditions was the last and 
long-standing problem of classical model atmosphere theory and it is finally 
solved with a high degree of sophistication. Application of these models leads 
to highly successful analyses of hot compact stars. Spectral properties from 
the extreme UV through the optical region are for the first time correctly 
reproduced by these models. The essential milestones for this development. 
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starting from the pioneering work of Auer & Mihalas [5] are: 

• Introduction of the Accelerated Lambda Iteration (ALI, or "operator split- 
ting" methods), based upon early work by Cannon [2] and Scharmer [3]. 
First ALI model atmospheres were constructed by Werner [8]. 

• Introduction of statistical approaches to treat the iron group elements in 
NLTE by Anderson [32,33]. 

• Computation of atomic data by Kurucz [38], by the Opacity Project [37] 
and subsequent improvements, and by the Iron Project [39]. 
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